1 Pre-requisitos e instalación de paquetes

Lista de paquetes espaciales en R https://cran.r-project.org/web/views/Spatial.html

install.packages(c("rgdal","rgeos", "maptools", "raster", "sp", "sf"))

Paquetes generales:

maptools: manipulación de datos geográficos
raster: Lectura, escritura, manipulación, análisis y modelamiento de datos espaciales en cuadrículas.
sp: clases y métodos de datos espaciales (Pebesma & Bivand, 2018).
rgdal: provee enlaces a la libreria geoespacial GDAL y acceso a las operaciones de proyección y transformación de la libreria PROJ4 (Bivand et al., 2017).
sf: combinación estandarizada de funciones en sp, rgeos y rgdal para la manipulación de datos espaciales vectorizados (Pebesma, 2017)

1.1 Organización y creación de proyectos en R

Guide to reproducible code (British Ecology Society)

Crear un nuevo projecto en R
Crear carpetas básicas (e.g., ./data, ./figs, ./functions, …)

└── my_awesome_project
    ├── src
    ├── output
    ├── data
    │   ├── raw
    │   └── processed
    ├── README.md
    ├── run_analyses.R 
    └── .gitignore

Fuente: https://www.r-bloggers.com/structuring-r-projects/

1.2 Ejercicio # 1:

Crear un projecto de R
Crear carpetas básicas para la práctica
Abrir el projecto
Crear un ‘R script’
Crear etiquetas como librerias y datos dentro del script (shortcut: Crtl+Shift+R)

2 Representación simple de datos espaciales

2.1 Tipos básicos de R

Números

1 or 2 or 4000 or 10e9...

Cadenas (or characters)

"Hola mundo"

Operadores lógicos

"TRUE" or "FALSE"

Vectores

c(1, 2) or c("uno", "dos")

Matrices

matrix(data = c(c(1,2), c("Bucaramanga", "Bogotá")), nrow = 2, ncol = 2)

##      [,1] [,2]         
## [1,] "1"  "Bucaramanga"
## [2,] "2"  "Bogotá"

Data frames

data.frame(Num = c(1,2), Nombre = c("Bucaramanga", "Bogotá"))

##   Num      Nombre
## 1   1 Bucaramanga
## 2   2      Bogotá

2.2 Puntos geográficos

x <- c(-73)
y <- c(7.13)

## Colombia continental
## Punta gallinas (lat= 12.46, lon = -71.67)
## Isla San José, Rio Negro Guainia Colombia (1.34, -66.7)
## Cabo Manglares (1.65, -79.01)
## Leticia (-4.23, -69.95)
plot(x, y, xlim = c(-79.01, -66.7), ylim = c(-4.23, 12.46))
# Bogotá
abline(h = 4.71, v = -74.07)

## Quibdo 5.69, -76.65 (max 851.7 mm)
## Medellin 6.24, -75.76 (226.7 mm)
## Bogota 4.71, -74.07 (119 mm)
## Bucaramanga 7.13, -73 (149 mm)
lon = c(-76.65, -75.76, -74.07, -73)
lat = c(5.69, 6.24, 4.71, 7.13)
ciudades = cbind(lon, lat)
precip = c(852, 227, 119, 149)

psize <- 1 + precip/500
plot(ciudades, cex=psize, pch=20, col='red', main='Precipitación', xlim = c(-79.01, -66.7), ylim = c(-4.23, 12.46))

https://rstudio.com/wp-content/uploads/2015/03/ggplot2-cheatsheet.pdf

library(tidyverse)
ciudades2 <- data.frame(x = lon, y = lat, precip = precip, nombre = c("Quibdo", "Medellín", "Bogotá", "Bucaramanga"))

ciudades2 %>% 
ggplot(aes(x = lon, y = lat)) +
  geom_point(size = 1 + precip/500) +
  xlim(c(-79.01, -66.7)) +
  ylim(c(-4.23, 12.46))

## Ejercicio

Busca en internet las coordenadas de 3 ciudades y la temperatura media anual
Crea un data frame con el nombre de la ciudad y las coordenadas
Crea una gráfica de los puntos con los límites deseados usando plot o ggplot

3 Conceptos espaciales

3.1 Vectores

3.1.1 Puntos geográficos

library(sp)
library(mapview)
library(tmap)

ciudades_geo <- ciudades2
coordinates(ciudades_geo) <- ~x+y

class(ciudades_geo)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

ciudades_geo

## class       : SpatialPointsDataFrame 
## features    : 4 
## extent      : -76.65, -73, 4.71, 7.13  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 2
## names       : precip, nombre 
## min values  :    119, Bogotá 
## max values  :    852, Quibdo

mapview(ciudades_geo, zcol = "nombre")

3.2 Proyecciones

Forma en que los datos espaciales son “aplanados” en un espacio 2D

“Alineación incorrecta de dos objetos espaciales = diferentes proyecciones”

Fuente: National Ecological Observatory Network (NEON)

3.3 Sistema de coordenadas

Fuente: Maning Sambale

proj4string(ciudades_geo) <- "+proj=longlat +datum=WGS84"
mapview(ciudades_geo, zcol = "nombre")

ciudades_geo

## class       : SpatialPointsDataFrame 
## features    : 4 
## extent      : -76.65, -73, 4.71, 7.13  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : precip, nombre 
## min values  :    119, Bogotá 
## max values  :    852, Quibdo

3.4 Extensión espacial

Coordenadas X y Y de las esquinas en el espacio geográfico

Fuente: National Ecological Observatory Network (NEON)

3.5 Polígonos

library(maptools)

## Checking rgeos availability: TRUE

data("wrld_simpl")

print(wrld_simpl)

## class       : SpatialPolygonsDataFrame 
## features    : 246 
## extent      : -180, 180, -90, 83.57027  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
## variables   : 11
## names       : FIPS, ISO2, ISO3,  UN,           NAME,    AREA,    POP2005, REGION, SUBREGION,      LON,     LAT 
## min values  :     ,   AD,  ABW,   4, Aaland Islands,       0,          0,      0,         0, -178.131, -80.446 
## max values  :   ZI,   ZW,  ZWE, 894,       Zimbabwe, 1638094, 1312978855,    150,       155,  179.219,   78.83

plot(wrld_simpl)

3.5.1 Filtración o subset (e.g., Colombia)

Col_inx <- which(wrld_simpl$NAME == "Colombia")
Col_poly <- wrld_simpl[Col_inx,]

plot(Col_poly)
plot(ciudades_geo, add = TRUE, col = "red", cex = 1 + ciudades_geo$precip/500, pch = 16)

Col_inx <- which(wrld_simpl$NAME == "Colombia")
Col_poly <- wrld_simpl[Col_inx,]

plot(Col_poly, axes = TRUE, col = "lightyellow")
plot(ciudades_geo, col = "orange", cex = 1 + ciudades_geo$precip/500, pch = 16, add = TRUE)

plot(wrld_simpl, xlim = c(-81,-67), ylim = c(-4,12), axes = TRUE, col = "lightyellow")
plot(ciudades_geo, col = "orange", cex = 1 + ciudades_geo$precip/500, pch = 20, add = TRUE)

plot(wrld_simpl, xlim = c(-81,-67), ylim = c(-4,12), axes = TRUE, col = "lightyellow")
plot(Col_poly, col = "orange", add = TRUE)
plot(ciudades_geo, col = "red", cex = 1 + ciudades_geo$precip/500, pch = 20, add = TRUE)

mapview(Col_poly, col.regions = "orange") + mapview(ciudades_geo, zcol = "nombre")

3.6 Cambiar la proyección de vectores

Lista de proyecciones disponibles

rgdal:::projInfo()

wrld_simpl2 <- spTransform(wrld_simpl, CRS("+proj=laea"))

par(mfrow = c(1,2))
plot(wrld_simpl)
plot(wrld_simpl2)

## Lambert Azimuthal Equal Area
ciudades_geo2 <- spTransform(ciudades_geo, CRS("+proj=laea"))

par(mfrow = c(1,2))
plot(ciudades_geo, col = "red", pch = 16, axes = TRUE)
plot(ciudades_geo2, col = "blue", pch = 16, axes = TRUE)

plot(Col_poly, col = "orange", axes = TRUE)
plot(ciudades_geo2, col = "red", cex = 1 + ciudades_geo$precip/500, pch = 20, add = TRUE)

Col_poly2 <- spTransform(Col_poly, CRS("+proj=laea"))

plot(Col_poly2, col = "orange", axes = TRUE)
plot(ciudades_geo2, col = "red", cex = 1 + ciudades_geo$precip/500, pch = 20, add = TRUE)

3.7 Ejercicio

Convierte el data frame del ejercicio anterior en un SpatialPointsDataFrame
Grafíca los puntos con diferentes proyecciones geográficas sobre un mapa de polígonos

3.8 Rasters

Datos espaciales en cuadrículas (“gridded data”)
Las cuadrículas representan areas de la tierra como pixeles

# <img src="http://www.neonscience.org/sites/default/files/images/dc-spatial-raster/raster_concept.png" style="width: 80%"/>

Fuente: National Ecological Observatory Network (NEON)

3.8.1 Resolución espacial

Tamaño de cada cuadrícula (usualmente en metros)

Fuente: National Ecological Observatory Network (NEON)

3.8.2 Datos climáticos

Variables bioclimáticas de la base de datos de Worldclim (Hijmans et al., 2004) http://worldclim.org/bioclim

BIO1 = Annual Mean Temperature
BIO5 = Max Temperature of Warmest Month
BIO6 = Min Temperature of Coldest Month
BIO7 = Temperature Annual Range (BIO5-BIO6)
BIO8 = Mean Temperature of Wettest Quarter
BIO9 = Mean Temperature of Driest Quarter
BIO12 = Annual Precipitation
BIO16 = Precipitation of Wettest Quarter
BIO17 = Precipitation of Driest Quarter
BIO18 = Precipitation of Warmest Quarter
BIO19 = Precipitation of Coldest Quarter

library(dismo)
# Descargar los datos directamente
bio_ly <- getData("worldclim", var="bio", res=10, path="./data/")

## Group of raster layers
class(bio_ly)

names(bio_ly)

tempRast <- bio_ly$bio1/10
precipRast <- bio_ly$bio12

tempRast
precipRast

http://worldclim.org/bioclim

# Desde una carpeta
library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

precipRast <- raster("./data/wc10/bio12.bil")
precipRast

## class      : RasterLayer 
## dimensions : 900, 2160, 1944000  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /home/se1212/Git_repos/UIS_course/Intro_Analisis_Espacial/data/wc10/bio12.bil 
## names      : bio12 
## values     : 0, 9916  (min, max)

3.8.2.1 Graficar raster

plot(precipRast)

## Usando la función spplot del paquete raster
spplot(precipRast)

plot(precipRast, xlim = c(-81,-67), ylim = c(-4,12), axes = TRUE)
plot(Col_poly, add = TRUE)

col_ras <- crop(precipRast, Col_poly)
plot(col_ras, axes = TRUE)
plot(Col_poly, add = TRUE)

spplot(col_ras, sp.layout = list("sp.lines", as(Col_poly, "SpatialLines"), col="white"))

3.8.3 Extracción de información de rasters

extract(precipRast, ciudades_geo)

## [1] 7399 1965  863 1204

ciudades_geo$precip_bioclim <- extract(precipRast, ciudades_geo)
ciudades_geo@data

##   precip      nombre precip_bioclim
## 1    852      Quibdo           7399
## 2    227    Medellín           1965
## 3    119      Bogotá            863
## 4    149 Bucaramanga           1204

mapview(col_ras) + mapview(ciudades_geo, zcol = "nombre")

3.8.4 Cambio de colores y “breaks”

library(RColorBrewer)

preciHist <- hist(log(precipRast), breaks=4)

preciHist

## $breaks
## [1]  0  2  4  6  8 10
## 
## $counts
## [1]   4732  23187 210982 334984   8996
## 
## $density
## [1] 0.004059148 0.019889995 0.180982053 0.287351964 0.007716841
## 
## $mids
## [1] 1 3 5 7 9
## 
## $xname
## [1] "v"
## 
## $equidist
## [1] TRUE
## 
## attr(,"class")
## [1] "histogram"

# Using the raster object with different breaks
plot(log(precipRast), 
     breaks = preciHist$breaks,
     col = brewer.pal(5, "Spectral"))

plot(log(col_ras), 
     breaks = preciHist$breaks,
     col = brewer.pal(5, "Spectral"))

plot(Col_poly, add = TRUE)

3.8.5 Ejercicio

Usando otra capa climática de bioclim (e.g., temperatura anual media):

Crear un mapa de con la capa climática usando la extensión espacial de Colombia
Extraer la información climatica de cada una de las ciudades usando el objeto ciudades_geo
Crear nuevos puntos de corte en el raster climático
Crear un mapa con los nuevos puntos de corte junto con los puntos geográficos de las ciudades

Objetos espaciales en R

Susy Echeverría-Londoño

1 Pre-requisitos e instalación de paquetes

1.1 Organización y creación de proyectos en R

1.2 Ejercicio # 1:

2 Representación simple de datos espaciales

2.1 Tipos básicos de R

2.2 Puntos geográficos

3 Conceptos espaciales

3.1 Vectores

3.1.1 Puntos geográficos

3.2 Proyecciones

3.3 Sistema de coordenadas

3.4 Extensión espacial

3.5 Polígonos

3.5.1 Filtración o subset (e.g., Colombia)

3.6 Cambiar la proyección de vectores

3.7 Ejercicio

3.8 Rasters

3.8.1 Resolución espacial

3.8.2 Datos climáticos

3.8.2.1 Graficar raster

3.8.3 Extracción de información de rasters

3.8.4 Cambio de colores y “breaks”

3.8.5 Ejercicio